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ABSTRACT 


iii 


M 87 has been observed in the 0.2-4 KeV X-ray band using the High 
Resolution Imager on the Einstein Observatory, and at 1.452 GHz using 
the Very Large Array. The radio map showed that the halo contained 
prominent asymmetries to the east and southwest. The X-ray map 
indicated similar asymmetries, but they were imbedded in the diffuse hot 
gas that surrounds the core out to a radius of several arcminutes. The 
hot X-ray emitting gas was assumed to be spherically symmetric and 
could, therefore, be subtracted from the image. The resultant image was 
asymmetric with major lobes to the east and southwest that coincide 
approximately with the asymmetries in the radio halo. The data 
indicates that inverse Compton emission is a plausible model for the X- 
rays coming from the asymmetric component. 
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Chapter I 
INTRODUCTION 


Purpose 

Inverse Compton emission is a basic physical process that involves 
the interaction of a low frequency photon with a relativistic electron. 
The electron gives up energy to the photon and consequently a higher 
frequency photon *s created. For ultrarelativi.stic electrons (Tucker 
1975) 



( 1 ) 


where Vi and Vf are the initial and final photon frequencies respec- 
tively, and y is the Lorentz factor. The value of Vf depends on the 
incident and scattering angles; and It has a maximum value of lJy 2 V£. 
Although the initial photons could be any frequency, in extragalactic 
astronomy one is frequently concerned with the situation where they are 
the 2.7K microwave background photons (vj = Hz), and the 
resultant photons are X-rays (vf « lO^-IO 1 ® Hz) when y is about 10 ^ to 
10 11 . 

The radio lobes of galaxies typically have magnetic fields of 
about 1 to 100 microgauss. When these same regions are emitting syn- 
chrotron radiation at radio frequencies (.01 to 100 GHz), then they are 
a source of ultrarelativistic electrons with a y of about 10^ to 10^. 
Therefore one might expect to see X-ray lobes of similar shape and 
extent as the radio lobes. That inverse Compton emission has not been 


2 




confirmed Is probably due to two factors. A combination of large mag- 
netic fields in the radio lobes, and low photon densities puts the 
expected inverse Compton flux beyond the limiting sensitivity of avai- 
lable X-ray detectors. 

The elliptical galaxy M87 has been shown to have an asymmetric 
radio halo (Kotanyi 1980; Andernach et al. 1979), with lobes extending 
several arc minutes to the east and southwest (see Figure 1). Fabri- 
cant, Lecar, and Gorenstein (1980) observed M87 in the 0. 7-3-0 keV X-ray 
band using the Imaging Proportional Counter (IPC) of the Einstein 
Observatory.) Their data did not show any asymmetries in the X-ray 
image. However, using the High Resolution Imager (HRI) on the Einstein 
Observatory, Schreier, Gorenstein, and Feigelson (1982) discovered 
asymmetries in the X-ray image of M87 to the east and southwest. The 
authors found a count excess above background of 0.04 HRI counts/second. 
Based upon the simple estimates of the magnetic field from radio maps, 
the expected inverse Compton contribution was about 0.02 HRI counts/ 
second, and therefore was considered a plausible emission mechanism. 

M87 is a good candidate for inverse Compton emission for a few 
reasons. First, the east and southwest asymmetries contained in the 
radio halo (Kotanyi 1980; Andernach et al. 1979) could be spatially 
correlated with the X-ray asymmetries found by Schreier, Gorenstein, and 
Feigelson (1982). Second, the magnetic field in the radio halo is low 
(about 4 microgauss), and the radio spectrum is steep with a = 1.2 
(Andernach et al. 1979). The purpose of this research, then, is to use 
Improved radio and X-ray data to: (1) create an X-ray map that shows 








the asymmetric component by removing the symmetric gaseous component; 
(2) calculate th© expected inverse Compton X-ray flux using the radio 
map and compare to (1); and (3) to determine whether or not the inverse 
Compton process is occurir.g. If so, then parameters such as the mag- 
netic field, and the slope and amplitude of the electron energy spectrum 
(Harris and Gtindlay ly/9) can be calculated more precisely. 

Inverse Compton Effect 

Scattering of the 2.7K background can be treated as Thomson scat- 
tering since hv* « m e c 2 , where v' is the Doppler shifted 2.7K radia- 
tion in the rest frame of the electron. The total ‘'verse Compton 
spectral power per unit volume is (Tucker 1975) 

j(V) « /P(V)N(Y)dY ergs cm" 3 sec” 1 Hz -1 , (2) 


where 


P(v) - Sirr^hc Jf(v/4Y 2 v i )n ph (u i )dV 1 (3) 

is the spectral power assuming Thomson scattering with a single elec- 
tron, N(y) is the electron energy distribution function, vj[ is the ini- 
tial photon frequency, n p h(vj) is the initial photon number density, 
and r 0 is the classical electron radius. Equation (3) was derived first 
in the rest frame of the electron and then transformed back to the 
laboratory frame. The function f(x) is given by 
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f(x) - x + 2x 2 ln(x) + x 2 - 2x 3 0<x<l , (4) 

and represents the Thomson spectrum for a single scattering. The vari- 
able x is the ratio of the emitted frequency to the maximum possible 
frequency 4‘if 2 Vj[. For synchrotron emitting radio lobes the electrons 
follow a power law of the form 


N(Y) - NY“ n . 


(5) 


Equation (2) can then be rewritten 

2 .4,1 r n + 4 n+ll ”1 

J<U>-*r ohclJ2 »« |(„ +3) 2 ( „ +lH „ +5) J 


( 6 ) 


v -(n-l)/2 j vj n_1 > /2 n ph (v 1 )dv 1 . 


The spectral index of the resultant radiation is (n-1)/2. Since the 
initial photons represent a blackbody, then (Rybicki and Lightman 1979) 


n P h ( V ’ 8lTV i = ” 3 (e h V i /kT - l)" 1 cm" 3 . (7) 

From equations (7) and (6), and equation 23.2.7 from Abramowitz and 
Stegun (1970), the final result is 

j(v) - 8<kT) 3 N(r o ir /he ) 2 ( kT/ hv) < n_1 > /2 b( n) (8) 

- 4.18xlO” 40 NT 3 (2.08xl0 10 T/v)* n " 1)/2 b(n) > 


pr 




where 
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b( n) - 2 


n+3 


?» +4n+ll 

(n+3) (n+l)(n+5) 


_ ,n+5v , , n+5 . 

r (“T”) 5("T") . 


(9) 


The values of b(n) are tabulated in Table 1. ?(x) is the Riemann zeta 

function 


5(x> “ “ k~* . 

k* 1 


( 10 ) 


The synchrotron radiation from the relativistic electrons within a 
radio lobe is given by (Tucker 1975) 


j(v)'- 1.70xI0' 21 Na(n)B (n+1)/2 (4.19xl0 6 /v) (n “ 1)/2 } OD 


whore 


(n-l)/2 


a(n) 


(3/tt) 


1/2 


r( 


3n- 1 
12 


) r ( 


3n+19. <£+5 


12 


) r (“f-) 


8(n+l) T (*2±i) 


( 12 ) 


The values of a(n) are tabulated in Table 1. These values were found to 
be incorrect (for n>4) in a number of sources, all of which referenced 
Ginzburg (1961J). The calculated values in Table 1 agree with those 
given by Tucker (1975). Appendix B contains the FORTRAN code used to 
calculate both a(n) and b(n). Equation (11) assumes the electron dis- 
tribution to be locally isotropic and follow the power law of equation 
(5). The electron spectral index, n, is related to the synchrotron 
spectral index, a, by the relation 
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Table 1 

The functions a(n) and b(n). The error is less than 1 x 10"^. 


a(n) b(n) 


0.5 

1.3635 x 10° 4 

3.0064 x 10 1 

0.6 

8.1296 x 10" 1 

3.0039 

0.7 

5.6718 

3.0239 

0.8 

4.2967 

3.0646 

0.9 

3-4275 

3.1253 

1.0 

2.8341 

3.2055 

1.1 

2.4068 

3.3052 

1.2 

2.0871 

3.4247 

U3 

1.8405 

3-5649 

1.4 

1.6460 

3.7267 

1.5 

1.4896 

3.9115 

1.6 

1.3620 

4.1209 

1.7 

1.2565 

4.3569 

1.8 

1.1685 

4.6218 

1.9 

1.0944 

4.9183 

2.0 

1.0317 _ 

5.2495 

2.1 

9.7823 x 10‘ 2 

5.6189 

2.2 

9.3255 

6. 0306 

2.3 

8.9341 

6.4892 

2.4 

8.5984 

6.9997 

2.5 

8.3107 

7.5682 

2.6 

8.0646 

8.2012 

2.7 

7.8551 

8.9062 

2.8 

7.6779 

9.6919 

2.9 

7.5297 

1.0568 x 10 1 

3-0 

7.4074 

1.1545 

3.1 

7.3088 

1.2635 

3.2 

7.2319 

1-3853 

3.3 

7.1751 

1.5215 

3.4 

7.1371 

1.6737 

3.5 

7.1167 

1.8442 

3.6 

7.1131 

2.0351 

3-7 

7.1256 

2.2492 

3.8 

7.1535 

2.4894 

3.9 

7.1966 

2.7592 

4.0 

7.2545 

3. 0623 

4.1 

7.3271 

3.4032 

4.2 

7.4134 

3.7870 

4.3 

7.5161 

4.2194 

4.4 

7.6326 

4.7069 

4.5 

7-7641 

5.2569 

4.6 

7.9109 

5.8781 

4.7 

8.0734 

6.5801 

*.8 

8.2519 

7.3741 

4.9 

8.4471 

8.2728 

5.0 

8.6596 

9.2909 


ft ‘ 

a ■ 

a 
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n - 2 ct + 1 . ( 13 ) 

From equations (8) and (11) it can be seen that the synchrotron and 
inverse Compton spectral indices are equal. 

For an optically thin source S v * j(v), where S v is the flux den- 
sity. Therefore 


S x /S r " J<v x >/j(v r > > 


(14) 


where S x and S r are the X-ray and radio flux densities. Using equations 
(8), (ID, and (13) this ratio becomes 


S /S - (v /v )“ 2.46 x10’ 19 t 3+ “b‘ ( “ +1) 


. (4 „96x10 3 )“ . 

a( n) 


(15) 


Because M87 has a low redshift (z = .004), corrections for the value of 
the temperature, T, would be negligible and have not been included. 


Morphology of M87 

M87 is a giant elliptical (El) galaxy, and the most massive member 
of the Virgo cluster. It has an apparent magnitude of about 10, and is 
15-7±0.7 Mpc distant (Mould, Aaronson, and Huchra 1980). The elliptic- 
ity varies from 0.76 to 3*72 at distances from the nucleus of 40" to 
32.5” respectively. An increase in ellipticity such as this is typical 
for elliptical galaxies. The position angle of the major axis is about 
150° and is fairly constant at all distances (Carter and Dixon 1978). 


M8T has proven bo be an interesting galaxy. The first extraga- 
iactic jet ever discovered was observed by Curtis (1918) on the nor- 
thwest side extending outward from the nucleus about 20". Moreover, M87 
is very rich with galactic clusters with current estimates at about 4000 
clusters (Harris and Smith 1976). At radio wavelengths, M87 is the 
third brightest source in the sky. The radio source identifications are 
Virgo A or 3C274. 

M87 was first detected in X-rays by Byram, Chubb, and Friedman 
(1966) using a rocket borne detector. Subsequent observations (Bradt et 
al. 1967, Adams et al. 1969, Bowyer et al. 1970, Lampton et al. 1971) 
confirmed the X-ray detection, however some confusion remained as to 
whether the X-ray emission was associated with M87 or with the center of 
the Virgo cluster. Using the Uhuru satellite, Kellogg et al. (1972) 
found the dominant X-ray emission to be extended by about 1° around, and 
associated with, M87. 

Models for the X-ray emission have been proposed using both ther- 
mal (Ginzburg 1967; Malina, Lampton, and Bowyer 1976; Lea et al. 1973) 
and non-thermal processes (Shklovskii 1967; Brecher and Burbidge 1972). 
Harris and Romanishin (1974) investigated several clusters of galaxies 
and found inverse Compton to be a viable mechanism for the generation of 
of the observed X-rays. It was Turland (1975) who first computed the 
expected inverse Compton X-ray flux from the radio halo. Using an 
equipartition magnetic field of 5 microgauss, he found the expected 
inverse Compton emission to be .02 percent of the total X-ray emission. 
Because thermal bremsstrahlung predicted the X-ray emission quite well, 
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Turland (1975) favored it as the emission process. The discovery of 
iron emission lines (Mushotzky et al. 1978, Lea et al. 1979) and 0 VIII 
emission lines (Canizares et al. 1979) confirmed the thermal bremsst- 
rahlung model. 

Whereas the early models of the M87 X-ray source were based upon 
low resolution data; Schreier, Gorenstein, and Feigelson (1982) used the 
Einstein HRI detector, which has a resolution of 3" full-width at half- 
maximum (FWHM) , to distinguish four regions: (1) a diffuse gas compo- 
nent; (2) a compact nuclear source; (3) an X-ray jet; and (4) an asym- 
metric component. The diffuse gas component is the large extended 
region that had been seen by earlier observers using low resolution 
data, and is the dominant or most luminous feature. The diameter of the 
nucleus is less than or approximately equal to 3", and may be associated 
with a central massive object. An accretion efficiency of < 10 -1 * onto a 
10^ M 0 black hole is required not to exceed the X-ray luminosity of lO 2 ^ 
ergs sec ■. The X-ray jet is very well correlated with the well defined 
jet observed at radio and optical wavelengths. The jet emission at 
X-ray wavelengths, as well as at longer wavelengths, is due to synchro- 
tron emission. The fourth component, the asymmetric component, is the 
one of interest for this research. This component is visible in Figure 
2 as deviations to the east and southwest from the otherwise spherical, 
concentric contours of the diffuse gas. The radial distribution of 
surface brightness is shown in Figure 3 (from Schreier, Gorenstein, and 
Feigelson 1982). Each set of points represents a 60° azimuthal bin with 
vertex at the nucleus. The enhancement in X-rays due to the jet can be 
seen around 10” from the nucleus in the 240° to 300° (plus signs) bin. 



X-ray map of the central ‘ 
The contour levels are 150, 82 
2.1, and 1.3 HRI counts. 
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Figure 2 




2' Gaussian smoothed with a FWHM of 20 
, 32.8, 16.4, 8.2, 6.3, 4.9, 3-3, 2.i 





The asymmetric component appears as a count excess from about 30” to 5' 
from the nucleus in the 60° to 120° (filled circles) and 180° to 240° 
(squares) bins. 


HRI BACKGROUND 


log r (orcsec) 


Figure 3 


Radial distribution of surface brightness about the nucleus in 60° 
azimuthal bins. The position angles defining the bins are 0°-60° (open 
circles), 60°-120° (filled circles), 120°-180° (crosses), l80°-240° 
(squares), 240°-300° (plus signs), and 300°-360° (triangles). The sub- 
tracted instrumental background is shown in the bottom left corner. The 
error bars at the top are the statistical uncertainties. The HRI point 
response function is shown by the dotted curve (from Schreier, Goren- 
stein, and Feigelson 1982). The surface brightness due to the diffuse 
gas component is plotted using stars. 
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Chapter II 
METHODS 

X-ray and Radio Data 

The radio and X-ray data for this research werss obtained in the 
form of a Flexible Image Transport System (FITS) tape. The FITS format 
is described by Wells, Greisen, and Harten (1981). 

The HRI was used to collect the X-ray data. The energy band in 
which the HRI responds is from 0.15 to 1.5 keV. The HRI itself has 4096 
by 4096 pixels, each subtending 0.5". The focusing optics of the Ein- 
stein telescopes limits the resolution to 3” (Giacconi et al. 1979). 
The X-ray data was made from two HRI exposures of M 87 . The first expo-* 
sure, which was used by Schreier, Gorenstein, and Feigelson (1982), was 
made on July 5, 1979 and lasted 75,600 seconds. The second exposure was 
made on December 9, 1980 with a duration of 51,200 seconds. The two 
images were coaligned and merged together. The data was reprocessed 
(Rev. 1.0) and only data corresponding to the best satellite aspect 
solutions (LOCK mode) were kept. The final effective exposure time of 
the merged image, i.e., the amount of time the instrument was actually 
collecting data, was 74,464.7 seconds. 

The X-ray image stored on the FITS tape was a 512 by 512 matrix of 
4" by 4" pixels. Rather than using the entire image, a 315 by 315 sub- 
set array contered on the nucleus was kept. This was done for two rea- 
sons. First, the subset array contains roughly one third the number of 
pixels in the full array. Therefore much less disk space was required 
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for storage; and the computation times were greatly decreased. Second, 
the subset array extends 10* from the nucleus. This is ample space 
since the radio lobes extend outward about 5'. Figure 2 is a contour 
plot of the central 12* of the X-ray data smoothed with a Caussian 
function with a FWHM of 20”. 

The radio data was collected from the National Radio Astronomy 
Observatory's (NRAO) Very Large Array (VLA), and represents the sum of 
two data sets. The first data set is a D configuration snapshot made by 
R. Perley in October 1981. The second data set is a hybrid C/A 
configuration synthesis obtained by E. Feigelson, E. Schreier, and D. 
Harris. For this configuration the telescopes on the north arm were 
more widely spaced than those on the east and west arms. Both data sets 
were observed at a wavelength cf 20 cm (1452.4 MHz). The self calibra- 
tion at high resolution of the second data set was not perfect. The 
data set contains a few residual stripes which show up as narrow linear 
contours. The wide spacing on the north arm and the related lack of 
north-south baselines, as well as the large dynamic range of radio 
intensities across the galaxy (roughly 1000 to 1 from the nucleus to the 
halo), are possible causes of the self-calibration difficulty. Despite 
these difficluties, the radio halo structure appears reliable. 

The radio data on the FITS tape was a 512 by 512 pixel array (each 
pixel was 4" on a side). Similar to the X-ray data, a 315 by 315 subset 
array centered on the nucleus was kept. Figure 4 is a contour plot of 
the central 12' of the data. 




Subtraction of the Diffuse Gas Component 
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j In order to calculate the X-ray contribution from the asymmetric 

I component, it was necessary to separate this component from the other 

| three components. The compact nuclear source is approximately 3" in 

j size or a little less than one pixel (4"). The X-ray jet extends for 

! 

about 20" or five pixels. Since the X-ray map is nearly centered on the 
nucleus (within 4" of the middle pixel), then ignoring the center pixels 
(out to a radius of five pixels) would be the simplest adequate way to 
remove the jet and nucleus. 

The diffuse gas component and the asymmetric component are not as 
easily separated because everywhere the asymmetric component occurs the 
diffuse gas component is also present. A simplification is possible if 
the diffuse gas is assumed to be spherically symmetric about the 
nucleus. From Figure 3 it can be seen that the slope of the different 
curves is approximately -1. Since it is a log-log plot, then a negative 
slope implies that the count density (HRI counts/square arcsecond) is 
inversely proportional to the distance from the nucleus. Therefore the 
; slope of the count density would be very steep close to the nucleus, and 

far away the slope would be near zero. The diffuse gas component can, 
with this 1 two assumptions above, be removed without need of any physical 
models for the type of emission involved. To do this, the X-ray map was 
divided into azimuthal pie slices and differentially spaced concentric 
circles, where two consecutive circles enclose an annulus (for example, 
the first and second circles define the first annulus). The first cir- 
cle had a radius of five pixels, therefore excluding the nucleus and jet 




fro* any bin. The basic idea is to make the bins small enough to accu- 
rately snaple the radial and azimuthal variations. Because azimuthal 
variations occur over large angles, a pie slice of 45° was adequate. 
Since this creates eight pie slices, they are hereafter referred to as 
octants. The circles were spaced according to the formula 

A ■ N/Rp j ( 16) 

where Rp is a constant, N is the circle number, and Afj+i is the 
increase in radius (which is the same as the spacing between circles) of 
the (N+1) th circle. The minimum value that An+; could be was set to 
1.5 pixels. The value cf Rp was chosen so that the spacing between 
circles, was small enough to accurately sample the steep gradient of the 
count density near the center, while at the same time avoiding needless 
oversampling further away where the gradient was small. Figure 5 shows 
the final binning of the X-ray map for Rp = 4 pixels (16"). 

Each annulus contains eight bins defined by the boundaries of the 
octants. The count density (here referring to HRI counts per pixel) was 
computed for each bin. Under the assumption that the diffuse gas com- 
ponent is spherically symmetric, the minimum count density within each 
annulus would represent that component without enhancement from the 
asymmetric component. Therefore, each annulus was assigned the lowest 
count density contained within it. Because each annulus was wider than 
one pixel, the minimum count density had to be associated with some 
specific radial distance. This distance was chosen to be the point that 
bisects the area of the annulus. Two arrays were, therefore, created; 
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OF POOR QUALITY 



Figure 5 


u and azimuthal binning of the X-ray map. Annuli 1 through 
12 have been left out for clarity. The filled circles indicate which bin 
had the lowest count density within that annulus. The underlying cont- 
our plot is the same as Figure 2. y 8 



one contained radial distances , and the second contained the minimum 
count densities associated with the radial distances. The minimum count 
densities are plotted on Figure 3 with stars. The points follow the 
lower edge of the data from Schrexer, Corenstein, and Feigelson (1982), 
which is the expected result if they are to represent the diffuse gas 
component. Therefore, using these points, and quadratic interpolation 
for intermediate pixels, the diffuse gas component can be subtracted 
from the X-ray image. The subtraction then leaves only the asymmetric 
component, which is shown in Figure 6. The FORTRAN program XSUB VSrOR 
(see Appendix A) performs the subtraction described above. It contains 
the variable ROTATE which rotates the pie slices from their zero posi- 
tion shown in Figure 5. To test the constancy of the subtraction, the 
pie slices were rotated 20° counter-clockwise. The subtracted data was 
found to deviate negligibly from the unrotated results. Therefore, the 
subtraction was invariant under this rotation. 
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Chapter III 
RESULTS 

The Radio Map 

The VLA data, used to create the radio map shown in Figure 4, is 
■not as good as the X-ray data from the HRI. Despite this, there are a 
few features that are probably reliable. The most prominent feature is 
the asymmetry in the radio halo. There is a distinct east lobe extend- 
ing about 4.5' from the nucleus, and a southwest lobe extending about 7' 
from the nucleus. These two lobes are seen in the data of Andernach et 
al. (1979) and the 610 MHz data of Kotanyi (1980). Their data is shown 
in Figure 1. There is also a smaller lobe extending 2" west of the 
nucleus, and directly opposite the east lobe. The smaller west lobe has 
approximately the same surface brightness as the larger east and south- 
west lobes. This feature does not show up in the data of Kotanyi 

(1980), but is somewhat visible in 4.8 GHz map of Andernach et al. 
(1979). A later 1465 MHz map by R. Perley is shown in Figure 7. This 
is a moderate resolution map (FWHM=62") which shows all of the features 
above as well as an elliptically shaped halo with a major axis about 15' 
long and elongated northeast to southwest. The southwest lobe is 
enveloped by the halo, whereas the east and west lobes appear distinct 
from it. 


The X-ray Map 

The asymmetric X-ray component shown in Figure 6 contains several 
interesting features. There are two lobes present. 




One lobe extends 





about 4* east from the nucleus, and has a FWHM of about 1' at a distance 
of 1' from the nucleus. The other lobe extends about 5' southwest from 
the nucleus. The FWHM of this lobe is about 0.7' at a distance of about 
2.5' from the nucleus. These lobes do not seem to extend any further 
than this, even at low surface brightness. This was indicated by the 
relatively steep gradient off the ends of the lobes on a more smoothed 
map. There is a "bridge" of emission extending about 2' south and 
southeast of the nucleus which connects the east and southwest lobes. 
Another feature seen is a large diffuse X-ray emitting region. It shows 
up as a sprinkling of low-level contours in the northwest and southeast 
quadrants, and extends about 15' from the northwest corner to the 
southeast corner. This large diffuse region could be the X-ray gas 
counterpart of the optical ellipticity. The position angle of the major 
axis is equivalent to that found by Carter and Dixon (1978) at optical 
frequencies. Figure 8 shows the above mentioned regions. Regions 1 and 
3 are the east and southwest lobes, respectively. Region 2 is the 
"bridge" of emission; and regions 4 and 5 together are the large diffuse 
region. 

In order to compare the X-ray asymmetric component and radio maps 
quantitatively, a conversion factor from HRI counts to Janskys is 
required. The values of % and a (the spectral index for the X-ray 
emission) need to be determined. From Stewart et al. (1984), the % 
contribution (assuming model C) from M87 was estimated to be 10 2 ^ cm“ 2 . 
However, the gas temperature is approximately 2 x 10?K which would com- 
pletely ionize elements up to about argon. Therefore the % contribu- 
tion from M87 would be negligible. From our galaxy % is generally 




accepted to be about 4 x 10 2 ® cm" 2 , and so this was the value used. The 
spectral index for inverse Compton radiation is equal to the spectral 
index for the synchrotron radiation of the relativistic electrons that 
are scattering the photons. Using the halo flux values at different 
radio frequencies from Andernach et al. ( 1979) > a was found to be about 
1 . 2 . 


The following equations convert HRI counts to flux and flux den- 
sity, 


NET/1000 - CTS * DTC * 1000 / EXP 


(17) 


F - CF p * 0.0 * 0.7924 * NET/1000 , (18) 


S “ CF S * 0.8 * 0.7924 * NET/1000, 


(19) 


where CTS is the total HRI counts, DTC is the dead time correction, EXP 
is the exposure time, CFp is the flux conversion factor, and CF§ is the 
flux density conversion factor. The conversion factors (for the given 
% and a) from the Einstein Observatory User’s Manual are, CFp = 1.47 x 
lO"^ in the 0.2 to 4 KeV band, and CF S = 1.61 x lO^ 1 at 1 KeV. Table 
2 gives the flux and flux densities for the five regions of Figure 8. 
The values given in Table 2 contain errors from a few different sources. 
One source of error is from the photon counting statistics. Because 
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Table 2 


X-ray fluxes and luminosities for the five regions of Figure 8. 
The area is in arcmin 2 j F x is in ergs cm -2 sec -1 ; L x is in ergs sec -1 ; 
and S x is in ergs cm -2 sec -1 Hz - '. 


Region 

Area 

F x (0.2-4 KeV) 

L x (0.2-4 KeV) 

S x (1 KeV) 

1 

6.12 

5.7x10 -12 

1.7X10 211 

6.2x10*3° 

2 

4.26 

2.5x10 -12 

7.3x10^° 

2.7x10 - 3° 

3 

5.88 

3.8x10 -12 

l.lxlO 1 * 1 

4.2x10*3° 

4 

23.4 

4.9x10 -12 

1.5x10^ 

5.4x10*3° 

5 

23-0 

4.0x10 -12 

1 .2x10^ 1 

4.4x10*3° 
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there are so many counts (on the order of a few thousand), a range of 
±1o represents less than ±2 percent of the total counts. This error is 
negligible compared to the other errors. An important source of error 
is in the conversion factors CFp and CFs. Given the uncertainties in a 
and Nfl, the conversion factors can vary by as much as a factor of two. 
In addition to this, the iumiridsifcy depends on tf»§ distance tg M87, and 
the value for this varies from 15*7±0.7 Mpc (Mould, Aaronson, and 
Huchra 1980) to 21.9±0.8 Mpc (Sandage, and Tammann 1976). 

Equation 15 relates the X-ray radio flux densities. It was modi- 
fied to use more familiar units. Putting in T = 2.7K and converting the 
units, the result is 


S../S. - ( v /E ) “ A .84xl0" 41 B~ ( “ +1) (5.54xl0 _5 ) a (20) 

a V n ) ) 


X r r X 


where S x is in ergs cm -2 sec -1 Hz -1 , S r is in Janskys, v P is in GHz, E x 
is in KeV, and B is in gauss. 
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Chapter IV 
DISCUSSION 


The results so far indicate that M87 contains an asymmetric radio 
halo which, because of the high surface brightness, is generally 
accepted to be emitting via synchrotron radiation. The X-ray image also 
contains an asymmetric component which is similar in size and shape to 
the radio halo. The emission mechanism for the X-rays could be linked 
to the ultrarelativistic electrons in the radio halo with an inverse 
Compton model. Alternatively, the X-ray emission could be due to a 
background cluster. 


, X-ray Cluster Emission 

Huchra and Brodie (1984) observed five galaxies that were within a 
few arcminutes of M87. Because their velocities are 25 times that of 
M87 indicates that they are background galaxies. The radial velocity 
distribution is about 300 km sec" 1 which indicates that they could be 
■embers of a poor cluster. The calculated distance to this possible 
cluster, using H 0 = 75 km sec" 1 Mpc" 1 , is about 343 Mpc. The positions 
of these galaxies with respect to M87 are shown in Figure 9. Two of the 
brighter galaxies, 87-1063 and 87-1169 (3 and 5 respectively on Figure 
9), coincide well with the southwest lobe. The other three appear to 
have no obvious connection with the asymmetric component. 

Huchra and Brodie (1984) suggest that the hypothesis of Schreier, 
Gorenstein, and Feigelson (1982) that inverse Compton emission is 


Figure 9 


Positions of the five background galaxies. Galaxy 87-869 is 1, 
87-522 is 2, 87-1063 is 3, 87-1303 is 4, and 87-1169 is 5. The thin 
contours are the asymmetric component (see Figure 6). The bold contours 
represent the poor cluster AWM 4 if it were behind M87 and 343 Mpc away. 



occuring in the halo requires reexamination in light of the discovery of 
this background cluster. In order to do this, an estimate of the clus- 
ter X-ray emission is required. Kriss, Cioffi, and Canizares (1983) 
made X-ray observations of 16 poor clusters. Unfortunately, the X-ray 
luminosity varied from < 1.07 x 10 1 * 2 to 1.02 x 10^ ergs sec” 1 , and so 
could not be used reliably. From the flux values for regions 1, 2, and 
3 (see Table 2), the expected X-ray luminosity for the asymmetric com- 
ponent, if M87 Here 343 Mpc away, is 1.7 x 10^ ergs sec” 1 . This is 
about 150 times brighter than the least luminous cluster, and about 1.5 
times brighter than the most luminous cluster. A convincing argument 
can be based on the shape of the asymmetric component. A "typical" 
X-ray cluster has a more or less spherical symmetry. The bold contours 
on Figure 9 are the X-ray map of a typical poor cluster, AWM 4, posi- 
tioned as if it were directly behind M87 and 343 Mpc away, which is the 
distance of the poor cluster discovered by Huchra and Brodie (1984). 
The figure shows that while the asymmetric component of M87 has about 
the right angular size, it is highly improbable that a background clus- 
ter would have such a unique shape and be centered around the nucleus of 
M87. From the above arguments it would be unlikely that the background 
cluster is contributing significantly to the observed X-ray asymmetry. 

The background cluster can also be shown to be unable to account 
for the halo emission at radio frequencies. Owen (1974) catalogued 
several hundred Abell clusters. Those clusters with a richness class 
R=0, and a distance class D=4, were used to obtain an average cluster 
flux density at 1400 MHz of 0.56±0.07 Janskys. This is negligible com- 
pared to the M87 radio halo flux of 80 Janskys at 1420 MHz. Therefore, 
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at radio wavelengths the background poor cluster could not be the source 
of the observed asymmetry. 

Inverse Compton Emission 

Because the radio data has a lower resolution than the X-ray data, 
the two maps could not be easily compared quantitatively using equation 
20. Qualitatively some interesting comparisons can be made. Figure 10 
shows the radio halo and X-ray asymmetric component maps plotted with 
the same scales. Both maps have major lobes to the east and southwest. 
The two east lobes have similar postion angles of about 85°. However 
the southwest lobe on the X-ray map has a position angle of 225°, while 
the southwest lobe on the radio map has a position angle of approxi- 
mately .210°. The small west lobe on the radio map appears to have no 
corresponding feature on the X-ray map. Similarly, the "bridge” between 
the X-ray east and southwest lobes does not appear on the radio map. 

The halo flux density could not be be accurately found from the 
radio data. The larger spatial features were not contained in the ori- 
ginal VLA database used to make the radio map shown in Figure 4. The 
result is that the data does not contain enough information to calculate 
the radio halo flux density. Therefore, to quantitatively check for 
inverse Compton emission, the radio halo flux densities of Andernach et 
al. (1979) were used. An equipartition magnetic field of 4 microgauss 
in the halo was assumed. This corresponds to a minimum energy for the 
relativistic particles and the magnetic field of about 10^® ergs. The 
relativistic particles are assumed to be entirely electrons (Andernach 
et al. 1979). The X-ray flux density used was the sum of the flux den- 
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slties from regions 1, 2, and 3 from Table 2. Therefore, using 
S x s 1.3 x 10” 2 ^ ergs cm” 2 sec” 1 Hz” 1 , E x = 1 KeV, B = 4 x 10"^ gauss, 
v r s 1,420 GHz, and a s 1.2, the expected radio flux density from equa- 
tion 20 is 128 Janskys, This is about 50 percent higher than the value 
of 80 Janskys given by Andernach et al. (1979). 

The calculation above involves considerable uncertainties, not the 
least of which is the magnetic field orientation. Equation 20 is the 
average over all angles. Dennison (1980) provides a very thorough ana- 
lysis of the problem of inverse Compton emission. He estimates the 
magnetic field in the halo to be about 1-5 microgauss, with the most 
probable value being about 2.5 microgauss. The expected inverse Compton 
contribution to the total X-ray luminosity in the 2-6 KeV band then is a 
little less than 1 percent. The total X-ray luminosity is about 
6.2 x 10 1 * 2 ergs sec” 1 (0.2-4 KeV). The asymmetric component (regions 1, 
2, and 3 from Table 2) has a luminosity of about 3-6 x 10* 11 ergs sec” 1 
(0.2«>4 KeV). Therefore these measurements indicate that inverse Compton 
accounts for about 6 percent of the total X-ray luminosity. 

Summary 

None of the data discussed so far provides conclusive proof for 
inverse Compton emission. However, calculated values are close enough 
to observed values that it is a very plausible mechanism. One important 
difficulty is that the percentage of inverse Compton luminosity con- 
tained in the total X-ray luminosity can vary from about 0.1 percent to 
more than 70 percent (Dennison 1980). The amount of the variation 
depends on the model, i.e., whether the X-ray emitting diffuse gas com- 
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ponent is isothermal or power law. There is also a strong dependence on 
the orientation of the magnetic field in the halo. A power law is the 
favored model for the diffuse gas (Stewart et al. 1984), and is consis- 
tent with inverse Compton contributing 6 percent of the total X-ray 
luminosity (see Figure 3» Dennison 1980). Therefore the inverse Compton 
intensity as seen by a particular observer can literally be turned on or 
off from region to region depending upon the orientation of the magnetic 
field, and the amount of random magnetic field that is present. This 
could, for instance, account for the spatial dissimilarities between the 
X-ray and radio map. These are not trivial problems to deal with. 
Their solution would appear to lie with better, higher resolution, radio 
data. Also polarization vectors for the radio halo would help to det- 
ermine. the direction of the magnetic field. 

At the upper limit, the X-ray cluster mentioned in a previous 
section could account for about 60 percent of the X-ray luminosity com- 
ing from the asymmetric component. This means that the cluster could be 
providing 4 percent of the total X-ray luminosity, and inverse Compton 
the remaining 2 percent. A 2 percent inverse Compton contribution is 
perhaps even more consistent with the consideratins brought out by Den- 
nison (1980). However, the unique shape of the asymmetric component 
would seem to preclude the X-ray cluster as a possibility. The cluster 
could have a luminosity closer to the lower limit, in which case the 
contribution would be insignificant. 

Region 2 (see Figure 8) is an interesting region. There is no 
corresponding radio feature visible on Figure 4. However, there is weak 
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radio emission from that region as can be seen on Figure 7. Region 2 
could be due to the background cluster, or a different emission process 
from inverse Compton. It is also possible that the X-ray emitting dif- 
fuse gas is centered southeast of the nucleus. 

in conclusion, the results of this work have reinforced the 
hypothesis of inverse Compton emission from the radio halo. They are 
consistent with the radio data of Andernach et al. (1979), but the 
inverse Compton X-ray luminosity is a little high according to Dennison 
(1980). The radio and X-ray asymmetric component maps show strong 
similarities, and the differences could be explainable as due to a 
change in the magnetic field's direction. 
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APPENDIX A 


FORTRAN Program bo Subtract the Diffuse Gas Component 


//V01XXXXX JOB 

// EXEC FVCG,PARM='0PTIMIZE(1)» 

/ /SOURCE. SYSIN DD * 

Cfttf«««***«««ft*»*«*ftft***«*«tt«ft*«tf***«ft*#**tt**##*#*»***tt*******»****##*«*ft 

C PROGRAM: XSUB VSFOR 
C 

C PURPOSE: To subtract the symmetrical (i.e., hot gas) component from 


C the HRI X-ray image of M-87. 

C 

C LOCAL VARIABLES: 

C ANN() Annuli boundaries 

C CTS() Counts within a given annulus and pie slice 

C DATA() Array to hold the data 

C DELPHI Size of pie slices in degrees 

C DENMINO .... Minimum count density for a pie slice in an annulus 

C FILIN Input file name 

C FILOUT Output file name 

C IDATA() ..... Array to hold position information 

C XYCEN, IXCEN . Center of symmetric component subtraction 

C JPHI Index number for pie slices 

C LBOX( ) Number of pixels within a given annulus 

C MAXANN Number of annuli 

C MAXREG Number of regions bounded by annuli 

C NBOX() Number of pixels within a given annulus and pie slice 

C NPHI Number of pie slices 

C NROWS,NCOLS . Number of rows and columns in the data 

C PHIMINO Pie slice containing DENMINO 

C PNAME Program name 

C R() The radius to point between annuli bisecting the area 

C RADDEG Factcr to convert radians to degrees 

C RADMAX Maximum distance from IYCEN, IXCEN for last annulus 

C RADMIN Minimum distance from center for first annulus 

C RDELMN Minimum allowed distance between annuli 

C RESID Hot gas component 

C RFACTR Factor used to determine the distance between annuli 

C ROTATE Rotation of pie slices in degrees 

C 


C Perry Wood, Eric Feigelson 
C Department of Astronomy, Penn State University 
C February, 1984 

C«**«tt«***«***tt«ft»«*****ft*»****ft*#**ftft#*********#*****ft***************** 

CHARACTER* 18 FILIN, FILOUT, PNAME 
INTEGER*2 NB0X(50,50), LB0X(50), IDATA(8) 

REAL*4 DATA(315,315), CTS(50,50), ANN(50), DENMIN(50) 

REAL*4 PHIMIN(50), R(50) 

C 

C . . . Initial parameters 
C 


anon 


C 


C 

C 

C 


DATA NROWS, NCOLS, IYCEN, IXCEN /315, 315* 158, 158/ 
DATA DELPHI, RADDEG, ROTATE /45., 57.29578, 20./ 

DATA FILIN, FILOUT /'XRAY.DAT01 1 , 'RESID.DAT02’/ 

DATA RADMIN, RFACTR, RQELMN /5., 4., 1.5/ 

DATA PNAME /'XSUB VSFOR'/ 


R1MAX s (NROWS - IYCEN) 
R2MAX = (NCOLS - IXCEN) 
RADMAX = AMIN 1 ( R 1 MAX , R2MAX ) 
NPHI = INT(360. / DELPHI) 

DO 14 I = 1 , 50 
DENMIN(I) s 1.0E20 
LBOX(I) « 0 
DO 16 J = 1, 50 
NBOX(I,J) = 0 
CTS(I,J) = 0. 

16 CONTINUE 
14 CONTINUE 


. . Read in the data 


DO 10 I = 1, NROWS 

READ ( 2 1 ) ( DATA ( I , J ) , J = 1, NCOLS) 
10 CONTINUE 

. READ(21)(IDATA(J), J = 1, 8) 

C 

C . . . Define annulus boundaries 
C 

DO 100 IRAD r 1, 5Q 
IF(IRAD .EQ. 1)THEN 
ANN(IRAD) = RADMIN 
ELSE 

DELRAD r IRAD / RFACTR 
IF(DELRAD .LT. RDELMN)THEN 
DELRAD = RDELMN 
END IF 

ANN(IRAD) = ANN( IRAD-1 ) + DELRAD 
END IF 

IF(ANN(IRAD) .GT. RADMAX) THEN 
MAXANN = IRAD - 1 
MAXREG = MAXANN - 1 
GO TO 110 
END IF 

100 CONTINUE 
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. . . Loop through the data array, and find in which radial and 
azimuthal bin each (IY,IX) pixel belong. 

110 DO 200 IY = 1, NROWS 
DELY = IY - IYCEN 
DO 210 IX = 1, NCOLS 
DELX = IX - IXCEN 
DIST = SQRT(DELX**2 + DELY**2) 

DO 215 IRAD s 1, MAXREG 




4 


ill 

IF(DIST .GE. ANN(IRAD) .AND. DIST .LT. ANN(IRAD+1 ))THEN 
JRAD = IRAD 

IF(DELY .LE. 0.0 .AND. DELX .GT. 0.0)THEN 
PHI = 270. - ATAN(DELY/DELX) * RADDEG + ROTATE 
ELSE IF(DELY .GT. 0.0 .AND. DELX .LT. 0.0)THEN 
PHI = 90. - ATAN(DELY/DELX) * RADDEG + ROTATE 
ELSE IF(DELY .GT. 0.0 .AND. DELX .GT. 0.0)THEN 
PHI = 270. - ATAN(DELY/DELX) * RADDEG + ROTATE 
ELSE IF(DELY .LE. 0.0 .AND. DELX .LT. 0.0)THEN 
PHI = 90. - ATAN(DELYZDELX) * RADDEG + ROTATE 
ELSE IF(DELY .GT. 0.0 .AND. DELX .EQ. 0.0)THEN 
PHI = 180. + ROTATE 

ELSE IF(DELY .LE. 0.0 .AND. DELX .EQ. 0.0)THEN 
PHI = 0.0 + ROTATE 
END IF 

IF(PH1 .GE. 360.) THEN 
PHI = PHI - 360. 

END IF 

JPHI = INT(PHI / DELPHI +1.) 

C 

C . . . Increment the count and area for bin (JRAD, JPHI) 

c 

CTS( JRAD, JPHI) = CTS( JRAD, JPHI) + DATA(IY,IX) 

NBOX( JRAD, JPHI) = NB0X( JRAD, JPHI) + 1 
LBOX(JRAD) = LBOX(JRAD) + 1 
END IF 

215 CONTINUE 
210 CONTINUE 
200 CONTINUE 
C 

C . . . Find and keep in array DENMIN(IRAD) the smallest count density 
C in each radial bin. 

C 

WRITE(6,6010)PNAME, RFACTR, ROTATE, FILIN, RADMIN, RADMAX, FILOUT, 
4 IYCEN, IXCEN 

WRITE(6,6020) 

DO 300 IRAD = 1, MAXREG 

R(IRAD) = SQRT( (ANN( IRAD)*#2 + ANN( IRAD+1 )#*2) / 2.) 

DO 310 I PHI = 1, NPHI 

DENS = CTS( IRAD, IPHI ) / NBOX( IRAD, IPHI ) 

IF(DENS .LT. DENMIN( IRAD) )THEN 
DENMIN(IRAD) = DENS 
PHIMIN(IRAD) = IPHI 
END IF 

310 CONTINUE 

WRITE{ 6 , 6000 ) IRAD , ANN(IRAD) , ANN( IRAD+1 ) , LBOX(IRAD), 

4 PHIMIN(IRAD), DENMIN( IRAD) , R(IRAD), R(IRAD)*i»., 

4 DENMIN(IRAD)/16. 

300 CONTINUE 
C 

C . . . Subtract the hot gas component. 

C 

DO MOO IY = 1, NROWS 

DELY = IY - IYCEN 
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420 


& 

& 


DO 410 IX = 1, NCOLS 
DELX s IX - IXCEN 
DIST = SQRT(DELX**2 + DELY*#2 ) 

RADM / IM ‘AMD. DIST .LT. R A DMA X) THEN 
IF(DIST .LT. R(2) )THEN 
I = 2 

ELSE IF(DIST .GE. R(MAXREG) )THEN 
I = MAXREG - 1 
ELSE 

J = MAXREG - 1 
DO 420 IRAD = 2, J 

IF I DI - ST IRAd' R(IRAD> ,AHD * DIST - LT * R(IRAD+1))THEN 

END IF 
CONTINUE 
END IF 

XXI = DIST - R( I— 1 ) 

XX2 = DIST - R( I) 

XX3 r DIST - R( 1+ 1 ) 

X1X2 = R( I- 1 ) - R( I ) 

X1X3 = R( 1-1 ) - R( 1+1 ) 

X2X3 = R(I) - R( 1+1 ) 

RESID = ( ( XX2 * XX3 * DENMIN(I-I)) / (X1X2 * X 1X3) ) - 
* XX3 * DENMIN(I)) / (X1X2 « X2X3)) + 

UXX1 * XX2 * DENMINd+1)) / (X1X3 * X2XR)) 
DATA(IX,IX) = DATA( IY, IX) - RESID 3>) 

r\ TP 


END IF 
410 CONTINUE 

400 CONS' 20 ’' 0 *™ 11 ' 1 *'' “ * '* BC0LS > 

WRITE ( 20 ) ( IDATA ( J ) , J = 1, 8) 

ST° p 

6000 F0RM4T( 2 X, I2, U.K.2, 14, 6X, FJ., ( J*, «.*, 

6010 F0RMAT(5X, 'Program: A18, T45, ' RFACTR = • F5.2 4X 

4 ^ROTATE = \F5.1, /, 5X, 'Input: A 18, T45, 'RADMIN = *, 

» 'ikuK’-.T V, 'ioM 1 : g; ;°) tput: ’’ fl ' 8 ’ Tfl5, 

6020 FORMAT ( ' IRAD ANN(IRAD) ANN( IRAD+1 ) LBOX 


£ 

i 


PHIMIN' 


END 


npiniTiM ^ , w .. ' ' ' rninirt 

98^?))’ 9X ’ R * ’OX, 'R'", 4X, 'CNTS/SQ ARCSEC 


W, 

/, ix, 


//DATA.FT21F001 DD UNIT=DISK,V0L=REF=MEN.U67410 PW6 LIB 
//^M=p4.U674IO.PH6.XRflY.Drai,DISp"(Oi:D,KEEP) ’ 
// DCBs ( RECFM=VBS , LRECL=200 , BLKSIZE=4004 ) 

//DATA.FT20F001 DD UNIT=DISK, VOL=REF=MEN.U674lO.FW6 LIB 

// JSN=MEN.U674l°. pw 6.RES I D.DAT02,DISP=(NEW,KEEpj, ’ ’ 

// DCB=(RECFM=VBS,LRECL=200,BLKSIZE=4004), 

// SPACE=(TRK, (20,10), RLSE) 


ooooooooooooo 
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APPENDIX B 


FORTRAN Program to Calculate a(n) and b(n) 


//V01XXXXX JOB 
// EXEC FVCG 
//SOURCE. SYSIN DD * 

PROGRAM: ABCALC VSFOR 


PURPOSE: Calculate the values of a(n) and b(n) from RADIATION 

PROCESSES IN ASTROPHYSICS, Tucker, equations 4-53 and 3-51. 

LOCAL VARIABLES: 


A a(n) 

B b(n) 

ERROR Desired accuracy of Riemann zeta function 


Perrv Wood 

lr 

Department of Astronomy, Penn State University 
July, 1984 

C******tf«*ft***»*«****««*******#»*#«*»«»****»**»#ft»*****»*«ft*****#»***** 

• LOGICAL* 1 RDATE(8) 

CHARACTER*8 WDATE 
CHARACTER* 15 PNAME 

REAL*8 A, N, C, B, N1, AB, BA, ERROR 
C 

DATA WDATE, PNAME, I VERS/* 07/3 1/84 ' * ‘ABCALC VSFOR', 1/ 

ERROR = 1.0D-6 
CALL DATE (RD ATE) 

WRITE( 6 , 604 ) PNAME , I VERS, WDATE, (RDATE(I), I r 1, 8), ERROR 
C = DSQRT(3.0D0) / (8.0 * DSQRT(3- 141592653589793DO)) 
WRITE(6,602) 

C 


C 


DO 10 N = 0.5, 5.1, 0.1 

A = 2**( (N-1 . )/2. ) * C * DGAMMA((3. * N - 1.) / 12.) 

4 * DGAMMA ( ( 3 . * N + 19 . ) / 12. ) * DGAMMA( (N +5.) / 4.) 

4 / ( (N + 1.) * DGAMMA((N + 7.) / 4.)) 

N1 = 0.5 * (N + 5.0) 

CALL RIEMAN(N1 , B, ERROR) 

B = B * 2.**(N+3.) * (N**2 + 4.*N + 11.) * DGAMMA(NI) 

& / ((N + 3-)**2 * (N + 1.) * (N + 5.)) 

AB = A / B 

BA = B / A 

WRITE(6 ,600)N, A, B, AB, BA 
10 CONTINUE 
STOP 


600 FORMAT ( 15X, F3-1, 4(3X, D12.5)) 

602 F0RMAT(16X, 'n', 8X, 'a(n)', 10X, 'b(n)', 9X, 'a(n)/b(n)', 

4 'b(n)/a(n) ' , /, 15X, ' — ', 7X, ' ', 1 1X, ' ' 

4 ' ' , 6X, ' •) 


6X, 

, 9X, 


1 1 " 




f 
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604 FORMAT(15X, 
& A8, 

END 


'Program: ', A15, 4X, 'Version: ', 12, 4X, 'Date: ', 
/, 15X, 'Run: 8A1, /, 15X, 'Error: < ', D 8 . 3 , //) 


C 

C 

C 

C 

C 

C 


This section calculates the Riemann zeta function. The sum is from 
zero to a cutoff point which is found by specifying the desired 
error. The error is less than or equal to the integration of the 
series from the cutoff to infinity. 

SUBROUTINE RIEMAN(S, ANS, ERROR) 

REAL*8 S, ANS, ERROR, CUTOFF, XERROR 
CUTOFF = 0.0D0 

20 CUTOFF = CUTOFF + 100. , k 

XERROR = ( CUTOFF** ( 1 .0-S) / (S - 1.0)) + CUT0FF**(-S) 

IF (XERROR .GT. ERROR) GO TO 20 
M = DINT (CUTOFF) 

ANS = 0.0D0 


DO 10 K s 1, M 

ANS = ANS + 1.0 / DFL0AT(K)**S 
10 CONTINUE 
RETURN 
END 


